A lattice gas model for active vesicle transport by molecular motors with opposite 

polarities 
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We introduce a multi-species lattice gas model for motor protein driven collective cargo transport 
on cellular filaments. We use this model to describe and analyze the collective motion of interacting 
vesicle cargoes being carried by oppositely directed molecular motors, moving on a single biofilament. 
Building on a totally asymmetric exclusion process (TASEP) to characterize the motion of the 
interacting cargoes, we allow for mass exchange with the environment, input and output at filament 
boundaries and focus on the role of interconversion rates and how they affect the directionality 
of the net cargo transport. We quantify the effect of the various different competing processes 
in terms of non-equilibrium phase diagrams. The interplay of interconversion rates, which allow 
for flux reversal and evaporation/deposition processes introduce qualitatively new features in the 
phase diagrams. We observe regimes of three-phase coexistence, the possibility of phase re-entrance 
and a significant flexibility in how the different phase boundaries shift in response to changes in 
control parameters. The moving steady state solutions of this model allows for different possibilities 
for the spatial distribution of cargo vesicles, ranging from homogeneous distribution of vesicles to 
polarized distributions, characterized by inhomogeneities or shocks. Current reversals due to internal 
regulation emerge naturally within the framework of this model. We believe this minimal model will 
clarify the understanding of many features of collective vesicle transport, apart from serving as the 
basis for building more exact quantitative models for vesicle transport relevant to various in-vivo 
situations. 



I. INTRODUCTION 

Biofilaments such as microtubules and actin play a cen- 
tral role in ensuring the mechanical integrity of eukary- 
otic cells [l| and intracellular transport. The entire net- 
work of actin and microtubule filaments serve as tracks 
for motor protein assisted transport of cargo vesicles and 
organelles from one specific location to another [l|, 
This transport process is active: multiple motor proteins 
bind to the filaments, and use ATP hydrolysis to convert 
the stored chemical energy into mechanical work 0, Q . 
In fact the heterogeneous distribution of organelles and 
vesicles and their regulation within the cell is achieved 
through this motor protein driven active transport pro- 
cess H, IU • The lack of detailed balance which under- 
lines the active motion of the molecular motors on fila- 
ments, renders the collective behavior of such a system 
markedly different from the one for an assembly of pas- 
sive molecules. 

Biofilament polarity provides a natural means to di- 
rect molecular motors, and depending on their molecular 
structure, different motor families displace to one of the 
two end of the filaments on consumption of ATP Q . For 
example, for microtubules (MT), dynein motors move to- 
ward the MT minus end while kinesin is a plus end di- 
rected motor The existence of motors moving in op- 
posite directions allows for bidirectional cargo transport. 
In-vitro experiments using fluorescence video microscopy 
with melanophore cellular extracts of cellular filaments ( 
actin and microtubules), motor proteins and vesicle car- 
goes [a [Z| has shown that under ATP rich conditions, 
individual cargo vesicles move bidirectionally on MT in 
the presence of oppositely moving motors such as dynein 



and kinesin [8j. 

A single cargo vesicle will be usually attached to several 
molecular motors If these motors have opposing 

polarity, the net cargo displacement along the filament 
will arise as a result of motor competition. Under these 
circumstances, interconversion or directional switching 
arising from motor competition and their kinetics to the 
filament track is required to ensure efficient cargo vesicle 
transport over long distances and avoid jamming. The 
origin of this competition, and the existence of regulatory 
mechanisms which determine the collective behavior of 
such motors remains an object of controversy [Tll | . While 
specific regulatory molecules have been proposed Q re- 
cent theoretical studies have shown that regulation can 
arise as a result of the non-linear dependence of the mo- 
tor velocity to the forces it is subject to [lS^ . Therefore, 
the traditional clear distinction between tug-of-war H, 
wher Even if insight has been gained in the molecular 
mechanisms underlying the transport of a single vesicle, 
much less is understood regarding bidirectional collective 
motion of vesicles [l4T j . Experimental evidence has shown 
that the spatial distribution of vesicles is sensitive to the 
cellular environment. Hence, their dispersion does not 
depend only on specific biochemical signals, which alter 
the motor binding properties, but can be controlled by 
modifying the polymerization degree of the embedding 
actin network. In experiments carried out with cellular 
extracts of fish melanophores, increasing the amount of 
caffeine, a chemical which alters motor binding kinetics, 
leads to a uniform dispersion of vesicle granules. On the 
contrary, adding latrunculin, which favors actin depoly- 
merization, the granules loose their ability to disperse 
and accumulate in the cell periphery [3, 6]. At high cargo 
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concentration, of relevance in in vivo situations, one can- 
not ignore inter-vesicle interactions. Studies carried out 
on mitochondrial distribution and movement in cultured 
neurons highlights the importance of inter-organelle in- 
teraction [15|. Moreover, transport is affected by cargo 
loading and removal at the filament end. For instance a 
sink mechanism has been proposed for transport of pig- 
ment granules in mammalian melanocytes, wherein gran- 
ules are captured by myosin V on reaching the filament 
end @, El]. 

These facts highlight the relevance of the interplay 
of different processes determining vesicle distribution. 
Coarse grained models which incorporate these compet- 
ing processes become useful means to rationalize their 
effects at large scales and describe the collective features 
resulting from bidirectionality. We will make use of such 
an approach to analyze quantitatively the emergence of 
spatial organization and regulation of vesicle distribu- 
tion induced by active motor transport. We will propose 
an effective model based on the totally asymmetric ex- 
clusion process (TASEP) [ljHij. TASEP exhibits non- 
equilibrium phase transitions between different macro- 
scopic density and current states. Variants of this model 
have been studied to model collective unidirectional mo- 
tion of molecular motors or other cargoes, whenever the 
detailed motor molecular interactions can be disregarded 
[21M241 ] . Despite its simplicity, this theoretical approach 
predicts inhomogeneous motor density profiles and shock 
development along the filament [2lj and provides a sim- 
ple reference framework to address generic issues related 
to collective motor transport [25[. TASEP has also been 
used to address the effective interaction with the filament 
environment, either because of diffusion or active trans- 
port through nearby filaments, through particle attach- 
ment/detachment, analogous to a Langmuir kinetic pro- 
cess (LK) . The phase diagram of TASEP-LK has shown 
the persistence of density shocks and the coexistence of 
different motor phases which characterize their collective 
transport behavior [26|, [2?J . 

Although, in general, a single vesicle cargo can be at- 
tached to a variable number of molecular motor, Ref. [l2[ 
has shown that the cumulative effect of competing mo- 
tors can effectively be described in terms of right- and 
left-moving species and that adding a resting species im- 
proves the comparison with experimental results on single 
moving vesicles [24| • For the sake of simplicity, we intro- 
duce a generalized two-species TASEP-LK model to de- 
scribe bidirectional collective cargo transport. Including 
a third, non-moving species, is straightforward and does 
not alter the qualitative conclusions we will discuss sub- 
sequently. The minimal model that we first put forward 
in Ref. 28] , accounts for the non-conservative and collec- 
tive transport of cargo vesicles on a biofilamcnt, which in- 
cludes the effects of excluded volume interaction between 
the cargo vesicles, evaporation-deposition processes of 
the vesicles in the bulk and boundary loading/off- loading 
of cargo vesicles at the boundaries. We will concentrate 
in a simple open geometry where vesicles enter from the 



extreme of the filament, as opposed to cyclic geome- 
tries, where for a single motor species a homogeneous 
motor distribution is expected. We will combine mean 
field analytical studies with Monte Carlo simulations to 
characterize the role of the boundaries in the structure 
and transport inside the filament. We will discuss when 
shocks develop and will determine the corresponding non- 
equilibrium phase diagrams, which show boundary and 
motor current reversals along with domain wall(shocks) 
localization along the filament. This simple theoretical 
framework allows us to identify the relevant physical pro- 
cesses which control the appearance and stability of the 
collective behavior of coexisting motors with opposed po- 
larity. Hence, despite its simplicity, the proposed frame- 
work will provide qualitatively understanding of the rele- 
vant features of bidirectional vesicle transport on micro- 
tubule networks. In Section [TT] we introduce the model 
we will study and in Section [IIII we derive the mean field 
(MF) continuum equations for the vesicle densities and 
describe also the numerical scheme we have implemented 
to check the MF theoretical predictions. Due to the open 
geometry under scrutiny, we analyze in detail the bound- 
ary conditions and their continuum limit. In Section ITVl 
we obtain and analyze the density and current profiles 
solving the MF equations obtained in the previous sec- 
tion. In Section [V] we construct the phase diagram after 
having obtained the corresponding phase boundaries. We 
follow it up with a discussion on the nature of the phase 
diagrams and conclude in Section IVT1 with a discussion of 
the main results and their implications for cargo vesicle 
transport. 



II. THE MINIMAL MODEL 

We represent the biofilament as a finite one- 
dimensional lattice of length L with N sites, labeled 
i = 0, . . . , N — I and the lattice spacing, e = L/N, is 
related the basic cargo vesicle displacement. The fila- 
ment ends are open and control the net cargo flux with 
the environment through specific kinetic rules. 

We characterize the microscopic state of the system 
at a given site i at time t in terms of the occupation 
numbers of the cargoes moving to the right, nf, and to 
the left, n~ . For microtubules they would correspond 
to the motion of cargo vesicles induced by kinesins and 
dynein motors, respectively. The cargo size prevents 
multiple occupancy, irrespective of the directionality of 
motion. Thus only one particle can occupy a given node 
at a given time so that the occupation numbers take 
values and 1. The filament voids occupation, n®, can 
be regarded as a useful auxiliary species, determined by 
the conservation law, nf + + ri" = I. 

For the sites i — 0, . . . , N — 2, right- moving parti- 
cles (+) can hop to the neighbouring site at i + 1 if it 
is empty with a rate kt- Correspondingly, for the sites 
i = 1, . . . , N — 1, left- moving cargoes (— ) can hop to the 
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by an interconversion rates, (fe-), corresponding to 
the rate at which a right (left) moving cargo becomes a 
lcft(right) moving one. These rates are considered uni- 
form along the filament. We display the relevant dynam- 
ical processes in Fig. Q] 




FIG. 1: (Color online) Schematic representation of all the 
dynamical processes for right (+) and left(-)-moving cargoes. 
The left and the right moving cargoes hop to the neighbouring 
site with rate 1 if the site is empty. The rates of entry, a± , and 
exit, (3±, of cargoes at the filament boundaries are displayed 
along with the adsorption, d±, and desorption, e±, rates for 
the left and right moving cargoes. At each node in the lattice, 
the interconversion rates from right (left) moving state to left 
(right) moving one are characterized by the rates k±. 



neighbouring site at i — 1 if the site is empty with a rate 
k t ; subsequently the inverse of this hopping rate will be 
taken as the unit of time. 

At the lattice ends, vesicle cargoes can enter and/or 
leave depending on their polarity. Therefore, at site i = 
0, right-moving particles (+) can enter the filament with 
rate a+, if the site is empty and left moving particles 
can leave it with a rate /3_ . At the opposite end, on site 
% — N — 1, (+)-particles leave the lattice with rate, /3+ 
and (-)-cargoes enter with rate a_ if the site is empty. 

Due to the finite processivity of the molecular motors, 
cargoes can detach from the filament, and diffusing vesi- 
cles can also approach and attach. The details of this 
interaction, however, depends on the properties of the 
medium. For example, a surrounding actin network al- 
lows for the jump of vesicles to the neighbouring tracks 
and the attachment is controlled by such active trans- 
port processes. We include vesicles hopping off and on 
the filament as kinetic attachment and detachment pro- 
cesses. For all sites in the bulk, i = 1,N — 2, (+) and 
(— ) moving cargoes can detach from the filament with 
site-independent rates e+ and e_ and attach at an unoc- 
cupied node with rates d+ and d- respectively. Finally, 
the interaction between motors of opposite polarity is 
accounted for by a local interconversion process between 
the two opposing species. We characterize such a process 



A. Relevant time scales and dynamic regimes 

Since there are different competing dynamical pro- 
cesses, qualitatively different dynamical regimes can be 
identified depending on the rate limiting mechanism. If 
all the rates exhibit the same scaling with particle num- 
ber, for a sufficiently long filament the dynamics corre- 
sponds to that prescribed by the attachment / detachment 
processes, and the behavior is characteristic of Lang- 
muir dynamics, except for a narrow region near the fil- 
ament ends [23]. This is the situation if processes are 
determined only by local interactions at the nodes be- 
cause the mean number of visited sites along the filament 
before a cargo detaches is of order l/e±. For attach- 
ment/detachment rates which scale extensively, a new 
dynamical regime is achieved where incoming fluxes com- 
pete with transport and desorption processes. An alter- 
native regime, where the environment exchange rates are 
dominant has also been explored in a related model for 
bidirectional motors (29j , which constitutes a generalized 
TASEP-LK scenario. We will concentrate in the regime 
for which the evaporation and adsorption rates scale 
inversely with the filament size; hence we introduce the 
rescaled rates E± — e±N and D± = d±N [23| . 
We will restrict ourselves to the scenario where intercon- 
version is the fastest process, relevant for bidirectional 
vesicular transport where opposite motors appear to en- 
gage in rapid succession In this limit, at each site 
i, the density of (+) species equilibrates with the density 
of (— ) particle and the two species become proportional 
to each other, as determined by the ratio of intercon- 
version rates. We will also see that in this regime the 
ratio of interconversion rates will determine the overall 
direction of the cargo flux. Moreover, the existence of 
two competing species breaks the particle/hole symme- 
try characteristic of one-species TASEP models, through 
the interconversion process. This lack of symmetry will 
affect significantly the non-equilibrium phase diagram of 
the system. 



B. Equations of motion 

Once we have established the dynamics of the two 
cargo species, we can write down the corresponding evo- 
lution equations for their occupation numbers. Assuming 
that the displacement rate is fast enough, we can treat 



4 



time as a continuous variable, leading to 
dn+ 



dt 



dn i 
dt 



n i-ii 1 ~ n i -n 4 ) - nf (1 -n i+1 - n 
kmf + kmf - e+nf + - nf - nf) 

n'i+xi 1 - n t - n f ) - - n \-i - 
fc 2 nf + kinf - e-nr + d-(l - nf -nf")(l) 



The terms on the r.h.s of Eq.([T]) describe gain and loss 
processes at a given node i arising from displacement, 
inter-conversion and deposition-evaporation events. As- 
sociated with displacement, we can identify the particle 
flux for the (+)- and (-)- particles at node i 



Jf = nf(l -n + i±1 -n j±1 ). 



(2) 



The boundary conditions are readily expressed in 
terms of the instantaneous vacancy number, nf At the 
left boundary site, i = 0, a vacancy enters, when either 
a (— ) particle leaves the left boundary or a (+) parti- 
cle hops to the neighbouring site in the bulk. Similarly 
a vacancy leaves the boundary site at the left boundary 
when either a (+) enters the boundary site, i = or a 
(— ) particle enters this site the next neighbouring site in 
the bulk, i = 1 Accordingly, 



dn° 
~dT 



nfn 1 + /3_n 



The terms on the right hand side of the equation above 
have a simple interpretation of gain and loss terms, with 
the first two terms contributing to the entry rate of va- 
cancies, while the last two terms correspond to the exit 
rate of vacancies at the boundary site. Similarly, the 
evolution of instantaneous vacancies at the right filament 
boundary reads 

dni 



N — l _ - , o + + 

— n N-i n N-2 + P+ n Ar-i — U-n N _ 1 — n N _ 2 n N _ 1 

(4) 

Averaging the occupation numbers over a time interval 
dt and over initial conditions, the corresponding equa- 
tions of motion for the expectation values are given by, 



d(nf) 



dt 



d(n~) 
dt 



= (nf-d 1 - nf -rii)}- (nf (1 - nf +1 - n l+1 )) 

- M n t) + M nj) -e+( nf) 
+ d+{l- n+ +1 - n- l+1 ) 

= K~+i(! -nf - rii)) - (nf (1 - nf_ x - nr_ x )) 



- M n i ) + k +( n t) -e-( rii ) 



(5) 



and accordingly, the mean void density at the ends of the 



filament can be expressed as 

^ = («)+/3_(n -}-a + (n0)-(nrn°) 



dt 



= ( n N-l n N-2) +P+i n N-l) - a ~( n N-l) 



( n N-2 n N-l) 



(6) 



The above set of equations are exact and, in principle, 
the steady states or the time evolution of the expectation 
values can be determined. However knowing the expec- 
tation value of any quantity requires the knowledge of all 
higher order correlations which makes it intractable ana- 
lytically [3(J. In specific cases, like TASEP, this moment 
hierarchy can be handled exactly using matrix method 
approach (l9j . Alternatively, if the number density fluc- 
tuations are smaller than mean density, a mean field ap- 
proach works remarkably well [20| 



III. MEAN FIELD (MF), CONTINUUM LIMIT 
AND MONTE CARLO (MC) SIMULATIONS 

A. Mean Field and Continuum limit 



We define the mean occupation densities as pf = (nf) 
and p~ = (nf ) , and we will introduce a mean field theory 
in which we factorize all the two-point correlators arising 
out of the different combinations of nf, nf as a product 



(3) of their averages, 



= nf)(n 



± 

i 

(nfnf +l ) = (nf)(nf +l ) 



± \ 



Pfpf+l 
Pfpf+l 



where the corresponding expression for the currents can 
be expressed as, 

Jf = (nf)((l-n + l±1 -n- l±1 )) = pf(l-pf ±1 -Pi ±i ) (7) 

If the filament length, L, is long compared with the cargo 
size, it is reasonable to keep L fixed and increase the num- 
ber of lattice sites, N — > oo, so that the lattice spacing 
e = jj — > [23[ and the position along the filament can 
be expressed as x = j^zf, with < x < 1. In this 
continuum limit, the mean densities can be written as 



(nf > 
( nf ±1 ) 



p+(x,t) 

p+{x,t)±e—p A 



(x,t) 



-—p +{x ,tl8) 



up to second order in e. Accordingly, the discrete mean 
field evolution equation, Eq. ([S]), with the corresponding 
boundary conditions, Eqs. (j6j), reduce to 



dp+ 
dt 



-P- 



eD+(l - p, 



- P- 

d 2 P . 

dx 2 



d 2 p + d 2 p- 



dx 2 dx 2 



P+ 



I - eE + p+ 

•(1-P+-P-) 

(9) 



5 



dp_ 

at 



—k-p- 



-P+ 



e£>_(l 
-P-)]- 



P- 
2 P- 



dx 2 



■(l-p+-p_) 



J_ = —K + J + . The total particle current along the fila- 
ment reduces to, 



d 2 p+ d 2 p- 



dx 2 



dx 2 



P- 



Rather than working with the two species densities, it 
will be useful to re-express Eq. © and Eq. (|10l) in terms 
of the total cargo density and the relative concentration 
of the different species, 



P = P^ 



P-> 



= 



P- 



(11) 



which yield 2p± = p ± <f>. Accordingly, we can rewrite 
the evolution equations to obtain 



dp 
at 



-£ = ^o(i-p) 



^[E+(P- 



00 

at 



2 [E+iP 



+ E-(p 
-E_{p-4>)] 



a 

e(D + 

a 



e 8x W 



(k- - k+)p 
(*++*-) 



K 



-p{x) 



(14) 



P+(x) = K+p-(x) 



J+(x)=p+(x) [l-W+p+{x)] 



(16) 



J(x) = J+(x) + J-(x) = (1 - K+) J+(x) 



(17) 



which explicitly shows that motor regulation allows for 
(lOcWrent reversal. For < K + < 1, the overall flux of 
vesicles is from left to right while for K + > 1, the overall 
flux is from right to left. Thus we see that for a fixed 
set of incoming and outgoing particle fluxes, a change 
in the motor expression can lead to total cargo flux re- 
versal. Different motor types exhibit different relative 
binding affinities Q. Accordingly, symmetric motor in- 
terconversion, K + = 1, when there is no net motor flux 
and cargo kinetics reduces to that of symmetric exclusion 
process [321 ]. This constitutes an exceptional situation, 
which we will not address in detail. 

The general expression for the local cargo density, 
Eq. (jT2j) , in the steady state reads 



(I-P)} (12) 
D-){l- P ) 
P)] (13) 



where we have introduced Dq = D + + £>_ , and have re- 
tained terms to leading order in e. These equations show 
that the dependence in terms of relative concentration 
is linear, and that non-linearities are only related to the 
global density occupation; the density dependence shows 
similarities with TASEP for a single species. 

We are interested in the steady state properties. Ac- 
cordingly, we set §2 = M = o. Since the interconversion 
rate is the fastest process, in the limit of small e, to lead- 
ing order the two densities become proportional to each 
other. 



where we have introduced the effective ratio of intercon- 
version rates K + = |p, which determines the amount 
of species imbalance at each node, and W+ = 1 + K + . 
From this relation, it is straightforward to obtain the cor- 
responding relation between the local occupancy of each 
species 



(15) 



which shows that K + , controls the local fraction of each 
species on the biofilament. For equal interconversion 
rates, the occupation is symmetric. 

The expression for current of (+) species reads in gen- 
eral 



Using Eq. (fT5|) it follows that the fluxes of left- and 
right-moving cargoes are proportional to each other, 



{2 P - i)d xP - (fi de + n ev ) P + n de = o (is) 

where, fide and fl ev are the total effective bulk deposition 
and evaporation rates respectively, which read 

Do _ E+ + E_ E+-E_ 

C ' ev ~ 2C + 2 



9. 



de 



(19) 



where, Co = ■ Eq. (fT5|) shows that the total den- 

sity evolves as an effective TASEP-LK model where the 
effective deposition and evaporation rates are modulated 
by the interconversion rate, K + . We can identify accord- 
ingly the effective binding rate constant, = Qde/Qev, 
which implies that in the absence of boundaries the 
cargo density relaxes to the Langmuir isotherm density, 
K lk 

P\k- 



l + K 



lk 



In the absence of interconversion, and for 



one moving species (i.e. K + = £L = D- = 0), K^. and 
Ph, reduce to those of single species TASEP, subject to 
Langmuir kinetics (23|. 

The particular case where fide = ^em -Ku, = 1, leads 
to symmetric behavior of the two motor species, where 
d _ , Ue + - E_) = A, and the overall cargo 



Co . 2Cn 1 2^+ ■ ; 

density follows the simple expression 



(2p-l) 



dp 
dx 



A 



= 0. 



(20) 



B. Boundary Conditions 

The presence of two species moving in opposite direc- 
tions alters the dynamics at the filament ends, whose 
behavior become relevant in the dynamical properties of 
ensembles of cargoes when the incoming and outgoing 
fluxes determine the particle distribution along the fil- 
ament. In the thermodynamic limit the entry and exit 
rates of the vacancies at site i = 0, read respectively 

Ren = p+(0)[l-p+(0)-p_(Q)]+/3_p_(0) (21) 
R ex - a+[l-p+(0)-p_(0)]+p_(0)[l-p + (0)-p_(0)] 
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At steady state R en = R ex . Using Eq. (IT5l) . the expres- 
sion for the density of (+) -movers at the left filament end 
is given by, 



P+(0) = 



M T [M 2 - Aa+W+ (1 - K+j\ ' 
2W+ (1 - K+) 



(22) 



where M = a+W+ + (3-K + + 1 — K + , which can also 
be derived in steady state from Eq.([6]). Similarly, at the 
other filament end, x = 1, 



P + [P 2 + 4a-W+ (1 - K+j\ 2 
2W+ (1 - K+) 



(23) 



where P = —a_W+ — f3 + + 1 — K + . The boundary 
densities do not depend on the attachment/detachment 
rates and coincide with the predictions for a bidirec- 
tional TASEP model decoupled from the environment, 
with D± = E± = HH. Unlike unidirectional TASEP- 
LK [23| , stochastic switching leads to boundary densities 
which do not depend linearly on the entry and exit par- 
ticle rates. 



C. Monte Carlo (MC) simulations 

In order to check the validity and limitations of the 
mean held approach introduced previously, we have car- 
ried out Monte Carlo (MC) simulations where we have 
implemented the dynamic processes described in Sec- 
tion HU A particle on the lattice and one of the different 
processes are selected at random for a MC move. A move 
for a particular process( i.e, translation or absorption- 
desorption) is accepted proportional to its rate. To 
achieve numerically the proper scaling between the rates 
associated with the different dynamical processes and to 
reach time scales in which cargoes can displace signifi- 
cantly along the filament, we take advantage of the fact 
that the interconversion dynamics (or stochastic switch- 
ing) between (+)- and (-)-motors is the fastest process. 
Hence, we enforce the time scale separation performing a 
large number of interconversion trial moves between each 
of the other two trial processes. We have checked that 
above 20 interconversion attempts pe 

We start from a random distribution of particles and 
let the system evolve before averaging. After an initial 
transient of around > 500^- swaps, where r stands for the 
rate of the slowest process (either deposition/evaporation 
or end filament particle fluxes), the system reaches its 
steady state. We then gather statistics of the relevant 
quantities averaging typically over 10 4 time swaps and 
collect information with a period > 10—. 



IV. DENSITY AND CURRENT PROFILES 

In order to establish the general non-equilibrium phase 
diagrams which characterize the collective motion of the 



two species of cargoes, we need to analyze the allowed 
density profiles in the regime where the dynamics along 
the filaments compete with the boundary fluxes. We will 
discuss separately the symmetric situation where the ad- 
sorption/desorption rates are equal, = 1, from the 
general case Kru ^ 1. The former, simpler case will help 
us to put forward the basic competing scenarios and to 
clarify the central role played by the incoming fluxes at 
the filament ends. 



We focus on the symmetric case to set the general 
structure of the collective motion of different species. Ac- 
cording to Eq.(l20p. the steady profiles read 



1 . .1 

p+ = oTtt - . J+ = P+0-- p+- P-) = 



2W 4 



(24) 



where use has been made of Eq. (|15p. The fluxes suggests 
that for the homogeneous profile, p = 1/2, the total cur- 
rent is maximized and thus corresponds to a maximal 
current phase, just like for TASEP system. We will refer 
subsequently to the homogeneous profile as p¥ . 

The inhomogeneous density profiles of + and — moving 
cargoes vary linearly and cannot satisfy the two bound- 
ary conditions at the same time. When the net particle 
flux is positive, K+ < 1, the overall density and individ- 
ual density for both (+) and (— ) species increases lin- 
early. We will distinguish between profiles which satisfy 
the density constraint at the left (right) filament end, 
i.e. low (high) density, p+ D (p+ D ). The corresponding 
phases will be named accordingly as the low density (LD) 
and high density (HD) phases. When the overall flux is 
from right to left, i.e when K + > 1, the solution which 
satisfies the density constraint on the right (left) bound- 
ary will correspond to the LD (HD) phase; the overall 
density profile grows from right to left. 

Without loss of generality, we will restrict ourselves to 
K + < 1 , because the behavior for K + > 1 can be derived 
under the transformation K+ O 1/K+, a+ <-> a_, /?+ 
/J_, E + <->■ E_, D + o Z?_ if we take into account that 
the overall flux reverses direction. Thus for K+ < 1 



P L + D {x) = p+(0) + B 3 



pl D {x) = p+(l)-B(l- : 



(25) 
(26) 



where B = -Dq/(1 — K + ). Here, p+(0) and p+(l) are 
the density of the (+) species at the left and the right 
boundaries, as prescribed by Eq. (f2"2"|) and Eq. (f2"3")l re- 
spectively. Although the rates of incoming and outgoing 
cargoes at the filament ends are natural control parame- 
ters, it is easier to classify the allowed coexistence regimes 
in terms of the densities at the filament ends. 

Let us consider a low density phase propagating from 
the left end of the filament and a high density phase 
propagating in the opposite direction while in the fila- 
ment center a maximal current phase can develop. If the 
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densities at the two filament ends satisfy p+(0) < 

and 2W~ — P+ W' we can identify the two position where 
the corresponding currents equal the maximum current. 
Such positions, x^ D ^ MC ^ and Xc° , respectively, 
determine the stability boundary of the corresponding 
phase. If the maximal current phase is involved, current 
continuity implies also a continuous change in particle 
density, similar to the situation in single species TASEP. 
However unlike TASEP, the LD and HD profiles are not 
spatially homogeneous. Using the current continuity con- 
dition, the crossover positions can be expressed in terms 
of the boundary densities as, 



x 



(LD/MC) 



1 



(HD/MC) 



2W+ 

(LD/MC) 



-p+(o) 



B 



(27) 



If both crossover positions satisfy < xi LD ^ MC ^ < 
x (hd/mc) ^ then the three regimes coexist in the 
filament, and the density profile is continuous and piece- 
wise linear for both species. We can then write down 
explicitly 



P+( x ) 



p+(0) + Bx 



(LD/MC) 



1 

2W+ 



p + (l)-B(l-x) 



, < x < x. 



(LD/MC) < K (HD/MC) 



(HD/MC) <x<1 



(28) 

If x^ D ^ IC ^ < 0, the low density phase is expelled from 
the filament. The inability of the maximal current phase 
to fulfill the boundary condition at the right filament end 
leads to a boundary layer. In the opposite case where 
x (hd/mc) ^ t w0 coexisting phases are the low den- 
sity and maximal current phases and a boundary layer 
develops at x = 1. Finally, if both crossover positions lie 
outside the filament length the maximal current phase 
occupies all the filament and boundary layers develop at 
both filament ends. In the limiting case where the two 

. .1 (LD/MC) (HD/MC) 

crossover positions coincide, x c — x c , the 

MC phase vanishes and particle densities move continu- 
ously from the low to the high density inhomogeneous 
phases and we have a two-phase coexistence regime. 
However for xi LD ^ MCS> > x ( ffD / MC ) f the intervening MC 
profile is expelled and ensuring current continuity for the 
LD and HD solution leads to a density discontinuity; 
hence the system develops a shock in the bulk, rather 
than at the filament ends. The shock position, from 
J+(x s ) LD = J+(x s ) HD , reads 



1 



X s — 



2B 

with a density jump 



H+ = p™(x s )-p^(x s ) 



1 

^LDi 



p+(0)-p+(l) 



(29) 



p+(l)-p+(0)-S (30) 
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FIG. 2: (Color online) Steady state spatial distribution of 
cargo density and current in the (LD-HD) phase coexistence 
regime, which exhibits a shock. Here a+ = 0.1, /3+ = 0.2, 
a- = 0.2, P- = 0.8 and D+ = 0.1, E+ = 0.2, D_ = 0.1, 
E- = 0.3 and K+ — 0.25. The dots represent the profile ob- 
tained by numerical simulations for (+)-moving cargoes while 
the squares are the profiles obtained by numerical simulation 
for ( — )-moving cargoes for a system size, iV = 1000. The 
lines represent the mean field solutions of the corresponding 
densities and current for (+)-moving cargoes which satisfy the 
appropriate boundary condition. 



provided < x s < 1. 



B. K lk >l 

Let us analyze the general situation when the symme- 
try between effective deposition and evaporation rates 
is broken. There exist two complementary strategies to 
determine the steady state density and current profiles. 
One can rewrite Eq. (|18|) in terms of p+, 



[2W+p, 



G P , 



B = 



(31) 



where G = D /(Cqp^.). Integrating Eq. (|3H and using 
the appropriate boundary condition at x = leads to 



2W + [ P+ (x)- P+ (0)} , 2W+B-G 



G 



Q2 



log 



Gp+(x)-B 



[G P+ (0)-B 



(32) 

and similarly the density profile consistent with the right 
boundary condition at x = 1 reads, 



2W+ 
G 

= x - 



[p + (x) - p+(l)] 

1 



2W+B - G 
CP 



log 



Gp+{x)-B 
G P+ {l)-B 
(33) 
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As usual, the boundary densities p+(0) and p+(l) are 
obtained from Eq. (JUJ) and Eq. ([25]) once the entry and 
exit rates for the two species are specified. Fig. [2] com- 
pares the MF density and current profiles and the corre- 
sponding MC simulation for a particular example, which 
shows density shocks in the bulk. 

It is however worthwhile to analyze the allowed st eady 
states in terms of the Lambert, W(x), function [27j j . 
This approach, which is discussed in Appendix A, is 
very useful in subsequently plotting and analyzing the 
resultant phases for our model. 

Depending on the boundary conditions, there exist three 
types of phases expressed in terms of the total density, 
p(x) = p+{x) + p-[x): a LD phase, Pl(x), with the 
density in the bulk is dictated by the left boundary 
condition; a HD phase, pr(x), and an analogous to the 
MC phase. As for the symmetric case, these phases may 
coexist, but the density profiles are no longer linear. 
Rather, the steady profiles approach monotonically the 
Langmuir isotherm density, py^. pl(x) never increases 
beyond 1/2, where the total current attains its maximal 
value. Similarly the solution, pr(x), always greater than 
1/2, approaches the Langmuir isotherm density, pri,, 
becoming independent of the boundary conditions from 
above (below) depending on the density of right moving 
motors, pr(1) > pru(pn(l) < puj |33j. Due to these 
properties, there are different scenarios depending on 
the imposed overall densities at the filament ends: 



If Pl(0) < 1/2,/9r(1) > 1/2, both solutions match 
the respective boundary conditions. HD and LD phases 
coexist and p(x) is obtained locating the position x w 
where the fluxes cross each other and the domain wall 
develops. If x w > l(x w < 0), no coexistence is allowed 
and the filament density is determined only by the left 
(right) filament end. 



When p L {0) < l/2,p R (l) < 1/2, the density profile 
matching the right boundary condition is unstable 
and the density approaches the extremal solution 
Wq(— Y*(x)), described in Appendix A, approaching 
1/2 at the right filament end. This phase, referred to 
as (M) phase [27j . is analogous to the maximal current 
(MC) phase described previously; the filament density 
is independent of the boundary value and the current is 
maximal. Unlike TASEP, the M phase has a spatially 
varying density profile , with the current approaching 
the extremal value only at the filament end. 
If pl(0) > 1/2, the bulk profile is already in the high 
density phase. However, for pr(1) < 1/2, the right 
solution is given as before by its extremal value and 
the entire bulk phase density ends in the M phase 
profile because matching the right boundary condi- 
tion is unstable. As a result, the profile at the right 
boundary approaches the extremal solution Wq(— Y*(x)). 



V. PHASE DIAGRAMS 

The analysis of the previous section provides a sys- 
tematic procedure to determine the allowed phases and 
derive the general non-equilibrium phase diagram for col- 
lective transport of competing cargo assemblies. Since 
the fluxes at the filament boundaries control the fluxes 
along the filament, the corresponding cargo rates are the 
natural parameters to describe the phase diagram. The 
absence of particle/hole symmetry leads to a non-linear 
relationship between particle fluxes and densities, intro- 
ducing significant differences in the phase diagrams with 
respect to those of single species TASEP. 
We will focus on phase diagram cuts in the (a + ,/3 + ) 
plane. We will discuss the main features of these phase 
diagrams and compare them with related lattice mod- 
els. Following the analysis of the previous section, we 
will first discuss the symmetric limit when the effective 
adsorption and desorption rates are equal. The insight 
gained will help us in addressing the general scenario. 



A. K lk = l 
1. The phase boundaries 

In the previous section we have identified the three 
relevant phases (LD,MC,HD), which can either cover the 
whole filament or coexist with each other. Hence, the 
non-equilibrium phase diagram can be built once we de- 
termine the limits of coexistence of these phases. To 
this end, we must determine the position of X ^ C LD ^ MC ^ ( 
^(hd/mc) an( ^ position of shock, x s and analyze 
when they are expelled from the filament. Hence we dis- 
cuss the coexistence onset between the different phases in 
terms of the boundary densities. The corresponding ex- 
pressions showing the explicit dependence on the bound- 
ary rates a+ and /?+ are shown in Appendix B. 

Phase coexistence line between (LD-HD ) and (LD-MC- 
HD): MC phase is expelled when x i LD ^ MC ^ > X \^- D I MC ^ 
Therefore, LD and HD phases will start to coexist with 
an MC phase when the density shock will vanish, i.e. 
^Jjld/mc) _ ^(hd/mc) a i on g fo e filament. In terms of 



the boundary densities, from Eq. (|2"T|) it follows 

p + (l)- P+ (0) = B, 



(34) 



from which the phase boundary in (oq_, /3+) plane is ob- 
tained by using Eqs. (|2"2"|) and ([2"3")l . which relate the 
boundary densities /Q+(l) and p+(0) with a+ and /3+. 
The obtained expression has two possible solutions de- 
pending on the parameters. Determining the physically 
plausible solution cannot always be resolved; we have 
then used Monte Carlo simulations to identify the ap- 
propriate boundary line. 

Phase coexistence line between (LD) and ( LD-HD ): 
When the shock position, x s , is close to the right filament 
end, then the profile in the filament is essentially LD 
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solution with an incipient HD profile. Accordingly, the 
LD-HD density shock starts at x s = 1. In terms of the 
boundary density, this condition reads 



P+(l) 



(35) 



Phase coexistence line between (HD) and ( LD-HD ) : 
Analogously to the preceding case, phase coexistence will 
start with a density shock at x s = 0, 



p+(l)+p+(0) 



1 



D 



(36) 



Phase coexistence line between (LD) and ( LD-MC ) 
and between (LD-MC) and (MC): When X [ LD/MC) is 
close to 1(0) the bulk profile corresponds to an incipient 
LD (MC) at the left (right) boundary and joining up with 
an emerging MC (LD) profile at the right (left) bound- 
ary. Therefore, the position of the coexistence line is 

1(0), 

(37) 
(38) 



determined by the limiting expression xi LD ^ MC ^ 
so that the equations of phase boundaries are, 

1 



Mo) 



2W+ 
1 

2W+ 



B 



for (LD)-( LD-MC ) and (LD-MC)-(MC), respectively. 
As shown in Appendix B, the coexistence curve is a 
straight line with fixed /3+. 

Phase coexistence curve between (HD) - (HD-MC) and 
between (HD-MC) - (MC): The coexistence starts when 
the MC phase is allowed at the right (left) filament end; 

n (HD/MC) =0{x{ HD/MC) = iy 



i.e. x> 



Mi) = 
Mi) = 



i 



2W 4 
1 

2Wl 



B 



(39) 
(40) 



for (HD)-(HD-MC) and (HD-MC)- (MC) coexistence on- 
set, respectively. 
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FIG. 3: (Color online) Phase space cut along the a+-/3+ plane 
for equal effective deposition and evaporation cargo rates ( 
A"^ = 1). (a) Without LK: When deposition and evaporation 
processes are not present then two or three phase-coexistence 
is absent and the phase diagram reduces to the case of two 
species TASEP with stochastic switching. Here /}_ =0.1 
, a_ = 0.3, K+ = 0.5. (b) LK is introduced, with Do = 
0.1 and the rest of parameters as before. There are clear 
topological changes in the phase diagram, with new regions 
of phase coexistence opening up. 
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2. Nature of the phase diagram 

Using the expressions for the coexistence curves, we 
can analyze the phase diagrams which emerge as a result 
of the transport of two competing species. 
The presence of two species moving in opposite direc- 
tions allow for overall current inversion as a function of 
the interconversion rate, a feature present in the absence 
of LK [3l[ . In Fig. [3] we compare the effect of LK in 
the collective cargo transport with respect to collective 
transport in the absence of LK. One can appreciate that 
interconversion prevents particle-hole symmetry, as ob- 
served above; moreover, LK is necessary to allow phase 
coexistence. Further, in the absence of LK the density 



FIG. 4: (Color online) Phase space cut along the a+-/3+ plane 
for equal effective deposition and evaporation cargo rates ( 
Auj = l)-(a) Without interconversion process (K+ — 0): 
Here (3- = , a_ = 0, D+ = 0.3, E+ = 0.3, D_ = 0, E- = 0. 
In this limit the phase diagram reduces to the case of single 
species TASEP-LK and the phase diagram corresponds to the 
ones reported in Ref. 27]. (b) K+ = 0.1, a_ = 0.1, /9_ = 0.3 
and Do = 0.3 and holding Ky^ = 1 . The phase diagram 
structure changes significantly from the previous case, with 
the LD region shrinking, accompanied by change in shape of 
(LD-HD) coexistence region. Also in this case the phase co- 
existence line between (LD-MC-HD) with (LD-HD) is not a 
straight line leading to the possibility of a re-entrant behaviour 
(see text). 
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and current profiles are homogeneous along the filament, 
unlike the ones corresponding to Fig. [2] 
In Fig. |4] we show the phase diagram in terms of the in- 
coming and outgoing rates of (+) species at the filament 
left and right ends, for a particular set of adsorption and 
desorption rates. As anticipated, by varying the parti- 
cle fluxes we can observe either density profiles corre- 
sponding to a single phase or coexistence of two or three 
different phases. Fig. [4ji shows collective transport in 
the absence of interconversion, when the model becomes 
equivalent to unidirectional TASEP-LK if a_ = /J_ = 0. 
Comparing it with Fig. 0b, one can appreciate that the 
breakdown of particle-hole symmetry induced by stochas- 
tic switching leads to a lack of symmetry in the phase 
diagram along the diagonal axis, a+ = /?+. This pecu- 
liarity leads to qualitatively new scenarios when motor 
species switch stochastically. 

For example, the phase boundary between (LD-HD)- 
(LD-MC-HD), shown in Fig. 0b and Fig. [5b are no longer 
straight lines, as is clear from Eq. (|B2[) . Hence, phase re- 
entrance is allowed; when moving along a linear path 
decreasing a + and increasing accordingly /3 + , it is pos- 
sible to pass from coexistence between a LD and HD 
profiles with a shock to a three phase coexistence (LD- 
MC-HD) and finally re-emerge into LD-HD phase coex- 
istence, exhibiting a re-entrant density shock (34j. As 
shown in Fig. 0b and Fig. [5b, topologically the phase 
boundary curve between (LD-HD)- (LD-MC-HD) region 
always ends at the vertex of another phase boundary 
curve. 
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FIG. 5: (Color online) Phase space cut along the a+-/3+ plane 
varying the total deposition cargo rates, Do = (D+ + D_), 
for fixed value of K lk = 2 and K+ = 0.2. (a) D = 0.02, (b) 
Do = 0.1, (c) D = 0.3, (d) Do = 0.9. E+ = E-, a_ = 0.2 
and /3_ = 0.8. For very low values of the total deposition rate 
as in (a), the topology of the phase diagram is similar to the 
one obtained in the absence of LK, where phase-coexistence 
is less favoured. 



B. Kyk >l 



We will follow the approach carried out in the previ- 
ous subsection and analyze first the coexistence boundary 
curves and analyze subsequently the main features of the 
corresponding phase diagrams. 



1. Phase boundaries 

Three phase coexistence is not allowed now due to the 
different nature of the maximal current phase. Therefore, 
the relevant curves reduce to three possibilities. 

Phase coexistence line between (LD) and (LD-HD): 
The incipient coexistence develops when the fluxes of the 
two phases match at the right filament end, leading to 



p L (l)+p R (l) = l 



(41) 



Relating the right boundary density with the particle in- 
put and output rates using Eq. (|23|) . one arrives at 



l-K, 



l-Pi(l) 



(42) 



in terms of pz,(l). Using Eqs. (|22p and (|A2I) one can 
derive numerically Pl(1) as a function of the input fluxes 



0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

p 1 l 1 l 




a 


i 


i 

M (<0 : 




\ 














t 


1 


HD 






i.i.i 



M W 



HD 



0.8 
H0.6 
0.4 
0.2 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

a + a. 



FIG. 6: (Color online) Phase space cut along the a+-/3+ plane 
varying for fixed value of K + — 0.2, E + = D_ = 0.05 
and a_ = 0.2 and /3_ = 0.8. (a) K lk = 1.2, (b)A" lk = 2, 
(c) Kn, — 6 and (d) Ky^ — 18. Increasing Ky^ decreases the 
low density region until it is eventually expelled; larger Ky^ 
implies higher bulk cargo density. 
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FIG. 7: (Color online) Phase space cut along the a+-/3+ plane 
varying the interconversion rate constant, K+, for fixed values 
of K lk = 2, E + = = 0.05 and «_ = 0.2 and /3_ = 
0.8. (a) K+ = 0.01, (b)K+ = 0.1, (c) K+ = 0.4 and (d) 
K+ = 0.8. The maximal current phase region is favored on 
increasing K+. The shape of the LD-HD coexistence region 
differs significantly from the shapes observed in the absence 
of interconversion 12711. 
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FIG. 8: (Color online) Phase space cut along the a+-/3+ plane 
varying the entry cargo rate, a_, of (— )-moving species for 
fixed values of K lk = 2, K+ = 0.2, /3_ = 0.8 and Do = 0.2. 
(a) q_ = 0.01, (b)a_ = 0.1, (c) a_ = 0.4 and (d) a_ = 
0.8. For high a_, system never attains current saturation for 
higher value of /3+ and in fact the M phase is absent, a feature 
missing in single species TASEP and TASEP with LK [13]. 
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FIG. 9: (Color online) Phase space cut along the q_-/3_ plane 
varying the interconversion rate constant , K+ , for fixed val- 
ues of K lk = 1.1429, a+ = 0.25, 0+ = 0.5 and D = 0.2. (a) 
K+ = 0.2, (b)K+ = 0.3, (c) K+ = 0.5 and (d) K+ = 0.7. For 
low K+ the system never attains current saturation and the 
M phase is expelled from the phase diagram 



After obtaining p L (0) using Eq. (J22J), Y L (x = 1) is 
evaluated. Using this, Eq. (|A2[) is solved numerically to 
obtain <tl(1) and hence Thus the entire phase 

boundary line in the a+ — /3+ plane can be constructed 
using these set of relations. 

Phase coexistence line between (HD) and (LD-HD): 
Phase coexistence starts when the density shock appears 
at the left filament end, x s = 0, leading to 



p L (0) + p R (0) = l 



(43) 



Using the relation of the left boundary density with the 
particle input and output rates through, Eq. (|22|) . one 
obtains the relation, 



ot+ 



T 2 - 2T(1 - K+ + (3_K+) 
2TW+ - 4(1 - K\) 



(44) 



where T — 2(1 — K + )p^(0) and we follow an analogous 
procedure to the previous case to relate pr(0) to the in- 
coming and outgoing cargo rates, a+ and j3 + . 

Onset of the M phase: The maximal current phase sets 
in when the maximum LK density, pr(1) = 1/2, develops 
at the right filament end. The corresponding value for 
the entry rate reads 



(3+ = a-W+ 



1 + K, 



(45) 



2. Nature of the phase diagram 



The topology of the phase diagram and even the nature 
of the phases are quite distinct from the symmetric case 
because three phase coexistence is not allowed. More- 
over, the phase coexistence lines separating LD and HD 
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have a different functional form and are not parallel any 
more. The nature of the M phase, although insensitive 
to the boundary rates, is no longer homogeneous and it 
is not characterized by a constant flux density. 
Figs. [5] and Fig. [6] illustrate how LD phase becomes less 
favored and it is eventually expelled on increasing the 
overall deposition, Dq, or K-^, leading to a topological 
change of the phase diagram. Such changes leave the 
M-HD coexistence curve virtually unaffected. Modify- 
ing instead the interconversion, K + , and (-)-cargo entry, 
a_, rates affect the position of the M-HD phase bound- 
ary as shown in Fig. [7] and Fig. [5] While always present 
in one species TASEP-LK [27|, the maximal M phase 
is destabilized on increasing a-, and beyond a critical 
value it is expelled from the phase diagram, modifying its 
topology as shown in Fig. [SJi. Interconversion enhances 
significantly the sensitivity of the distance between the 
(LD)-(LD-HD) and (HD)-(LD-HD) enhancing the corre- 
sponding coexistence regions, as displayed in Figs. [7^, and 
[7b. Fig. [7J shows how distance between these two coexis- 
tence lines varies significantly and vanishes for very small 
values of a+ and (3+ . 



3. Phase diagram in terms of incoming and 
outgoing rate of ( — ) species 

The expressions derived above describing the phase 
boundaries can be used to describe the phase diagram in 
terms of different control parameters. For example, if we 
use the entry rates of the (— ) motors, we arrive at anal- 
ogous phase diagrams, as shown, for example, in Fig. [5] 
Although similar to the phase diagrams in terms of the 
entry rates of the (+) species, the phase boundaries are 
rotated. For example in Fig. [H] the bottom right corner is 
HD phase, the bottom left is M phase, the middle region 
is essentially a LD-HD coexistence region and the top left 
corner is LD phase. Due to the symmetry between the 
two motor species, the phase diagram in a+, /?+ is identi- 
cal to the corresponding one in the a_,/3_ plane for the 
transformation K + — > K- = fc_/fc + ), E + — > E_, 

D + ->• D-, E- -)■ E + and £>_ D + , although with 
a reversed overall motor flux density. It is worth noting 
that the diagrams in the a+ , /3+ plane for K + < 1 are 
related by mirror inversion and a 90° rotation, with the 
corresponding diagrams in the a_,/3_ plane for K + < 1 
( and equivalently with plane for K + > 1). This 

symmetry is appreciated by comparing the phase plane 
plots of Figs. II and M 



VI. CONCLUSIONS 

We have presented a minimal model for bidirectional 
cargo transport on biofilaments to address the basic fea- 
tures associated to their collective dynamics. We have 
identified the relevant competing rates, the central role 
played by motor regulation, and have described the ex- 



pected scenarios depending on their relative magnitudes. 
We have shown that despite their opposed directionality, 
the collective transport may favor one particular global 
direction, in agreement with experimental observations. 
We have analyzed how current reversal emerges and its 
impact in cargo density accumulation at either end of 
a polar filament. We have focused on the regime where 
interconversion is the rate limiting process, a regime con- 
sistent with experimental observation of vesicle trans- 
port [ll|. The rate of interconversion may affect the 
phase diagram shape, as mentioned in Appendix C, and 
the framework described can be easily adapted to ana- 
lyze such regimes. 

This 'fast' interconversion regime, where the interconver- 
sion is the rate limiting process is simpler in nature and 
it has served to pinpoint the main features of collective 
transport of opposed polarity species. The derived non 
equilibrium phase diagrams rationalize the interplay be- 
tween motor regulation and motor exchange and they 
exhibit interesting and rich behavior and shows the rele- 
vant role played by motor regulation when several motor 
species interact to promote active cargo transport. In 
particular, we have described how shock re-entrance is 
allowed, and also how interconversion destabilizes and 
prevents certain phases for appropriate parameters. The 
interplay of the different processes determining collec- 
tive cargo motion provides a large degree of flexibility in 
phases and structures which cargoes can develop as com- 
pared to cargoes pulled by one molecular motor species. 
For example, both the maximal and low density phases 
can be expelled, and the regimes of LD and HD pro- 
files become more common as compared to one species 
TASEP. Moreover, the regulation of the two opposed 
species can lead to flux reversal without modifying the 
underlying phase diagram. 

Despite the simplicity of the model, we can relate the 
relevant model coefficients K+, ifvu and Dq to exper- 
imentally relevant situations. In-vitro experiments with 
cellular extracts of MT and molecular motors, motor reg- 
ulation (and thus effectively K + ) is achieved controlling 
dynein cofactors pl50/G/ued and dynamitin concentra- 
tion, which are subunits of the dynactin and regulate 
the relative binding affinities of dyneins to the MT Q. 
Experimentally it is known that increasing the cofactor 
concentration leads cargo motion reversal [3|, [6| . K-^ and 
Dq can be related to the motor interaction with the fil- 
ament and can be controlled by modifying either motor 
processivity d, 0] , filament concentration or the attach- 
ment properties of the cargo to other filaments. 

In order to address the competition between different 
mechanisms and how they affect the collective properties 
of cargo transport, we have analyzed how different phases 
develop along a filament as the attachment/detachment, 
interconversion and boundary input (output) rates vary. 
In particular, we have shown how one can move from a 
state of (LD-HD) phase coexistence, characterized by ex- 
istence of density shocks, to states with smooth density 
profile varying either the interconversion rate, K + , or the 
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effective desorption rate, D* . Therefore, the interplay 
between two different regulatory mechanisms, the detach- 
ment/attachment processes( e.g. due to the presence of 
an embedding actin network) and interconversion due to 
motor regulation, can give rise to a rich variety of inho- 
mogeneous collective cargo distributions. Understanding 
the impact of such scenarios in realistic biological systems 
remains an open challenge. 
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The relevant branch of the function is selected through 
the boundary conditions. Specifically, if the density at 
the left boundary lies in the range < p(0) < 1/2, then 

a L {x) =W- 1 {-Y L {x)) <0 (A6) 

while the solution which satisfies the right boundary con- 
dition at x = 1 fulfills 

(W (Y R (x))>0 , PR (l)> Plk 
a R (x) = { (A7) 
{ W (-Y R {x)) < ,p R (l) <p lk 

where Y^{x) and Y R (x) are the functions Y(x) introduced 
in Eq. (|A3[) obtained by setting xq = and xq = 1 re- 
spectively. 



Appendix A: Analysis of steady states using 
Lambert function 

It is useful to analyze the allowed steady states intro- 
ducing the auxiliary function a(x) (27[, 



a(x) 



K 



lk 



-[2p{x)-l]-l 



(Al) 



where p(x) = p+(x) + p~{x) refers to the total density. 
Expressed in terms of <j(x), Eq. (|18|) can be integrated, 
leading to 



|cr(x)| exp[cr(x)] = Y(x) 



(A2) 



where, 
Y(x) = \a(x )\ exp 



(K + 1) 

Qev—r; - Xq) + Cr(x ) 

^lk - 1 



(A3) 

Here, xo, refers to the position of the left (right bound- 
ary and cr(xo) is the corresponding value of the function 
at that boundary. 

Eq. (|A2[) has a explicit solution in terms of the Lam- 
bert, W(x), function [27j]. 



a{x) 



W(Y(x)) ,ct(x)>0 
W(-Y(x)) ,a(x) < 



(A4) 



The Lambert function has two real branches, Wq(Y(x)) 
and W-i(Y(x)), in terms of which we can write 



a{x) 



-2K U 

W-^-Y) :p€[0,±] )f7 g^Jk -1] 

Ik 



Wo(-Y) ,,£[-1,0] 



W (Y) :p e[^t_,l] ,^[0,^] 



(A5) 



Appendix B: Coexistence boundary lines 

In this appendix we provide the explicit expressions for 
the coexistence boundary lines in terms of the entry and 
exit cargo rates j3± and a± for symmetric TASEP-LK, 
where = 1. In this way, we complement the infor- 
mation discussed in Section [V A 11 providing explicit ex- 
pressions for the phase boundary lines used to construct 
the phase diagrams shown in Figf?] and FigOJ In order 
to proceed, it is convenient to express the boundary lines 
in terms of the following parameters: 

Do 



V 
Bi 

C 

Q 

Zi 

z 2 
z 3 
z 4 
z* 



1-K+ 
4a_(l -K%) 
1-K+ + /3-K+ 

V 

Mi-Kl) = — 

2{d + po)(l-K%) 
1 



= 2( 



1 



d-po)(l-K±) = -Zi + 



2W+ 



2(^ r +d-p )(l-K 2 + ) 



1 

~2W + 



(Bl) 

Here po = p+(0), which is given by Eq. (j2"2")l . Since p+(0) 
depends only on entry rate of (+) particle and exit rate 
of (— ) particles, Z\, Z 2 , Z%, Z4 and Z5 do not depend on 

P+. 

1. Phase boundary between (LD-MC-HD ) and (LD- 
HD): We solve Eq. using the expression for the 
boundary densities in terms of the incoming and 
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outgoing rates, Eqs. (1221) and ([23]) . The relevant 
rates are f3 + and a+. The phase boundary line in 
plane is given by, 

P+ = (l-K + -a-W+) + ^^, (B2) 

Note that the right hand side of Eq. (|B2[) is an 
explicit function of a + . Unlike the symmetric case, 
-RTu, = 1, the entry and exit rate, /3+ and a+, arc 
not proportional to each other. 

2. Phase boundary between (LD-HD ) and (LD): In 
this case, the emergence of coexistence is deter- 
mined by Eq. (I35|) using the explicit expressions 
for the boundary densities (Eqs. (|2"2")l and ([2"3l) h 

/3 + = (l-K+-a-W+) + ?-^-, (B3) 

3. Phase boundary between (LD-HD) and (HD): The 
coexistence curve is in this case encoded in Eq. (|36|) 
and on the expression of boundary densities as be- 
fore. The explicit expression is given by, 

i 8 + = (l-K + -a_W + ) + ^—?l, (B4) 



6. Phase boundary between (HD) with (HD-MC): 
This boundary is characterized by Eq. (|3"5)) and de- 
pends only on the cargo density on its right filament 
end. Using Eq. (|2"3"j) 

f3 + = (l-K + -a_W + ) + ^^- (B8) 

This phase boundary is a straight line parallel to 
the a+ axis. 

7. Phase boundary between (MC) with (HD-MC): Fi- 
nally the boundary line between these two phases 
obeys Eq. (|40|) and and depends only on the cargo 
density on its right filament end. Using Eq. (|23l) . 
we can explicitly write the expression for the phase 
boundary as a function of /3 + alone. 

f3 + = (l-K + -a_W + ) + ^^- (B9) 

This phase boundary is also a straight line parallel 
to the a+ axis. 



Appendix C: Deviations from the fast 
interconversion limit 



4. Phase boundary between (LD) and (LD-MC): Now 
the phase boundary is determined by Eq. (|37| . 
which depends only on the cargo density on its left 
filament end. Using Eq. (|22|) and introducing the 
parameter, 



c = ®(— 

1 2 V 2W- 



-d) 



the equation of the phase boundary reads 

C\ — 2B\C\ 
a+ = 2dW+ - Q 



(B5) 



(B6) 



This phase boundary is a straight line parallel to 
the /3+ axis. 

5. Phase boundary between (MC) and (LD-MC): The 
phase boundary is characterized by Eq. ()38|) , which 
also depends only on the cargo density on its left 
filament end, as in the previous case. The equation 
of the phase boundary is, 



C 2 - 2B X C 



(B7) 



2CW+ - Q 

This is also a straight line parallel to the /3+ axis. 



Although we have focused in the limit where effective 
motor interconversion is the fastest process, we have car- 
ried out some Monte- Carlo studies to analyze the impact 
of such an assumption on collective motor dynamics. We 
compare the interconversion with the translation rate, k t 
through the parameter q = k + /k t , which is infinite in 
the fast interconversion limit. For example, for a given 
set of parameters where the density profile is homoge- 
neous , as we decrease q we move to a regime where the 
density profile is inhomogeneous and where its details de- 
pend separately on the two interconversion rates, k + and 
In this regime, the simplifying assumption of pro- 
portionality between the density of right and left moving 
cargoes no longer holds and a separate treatment for the 
two species will be required. Interestingly, decreasing fur- 
ther q, we find typically a different homogeneous profile 
where the two cargo densities ( the left and right mov- 
ing cargo densities) become proportional to each other 
again, where the profiles depend only on interconversion 
rates through K+ . Fig. [10] shows a particular example 
where we start from an LD profile in the fast intercon- 
version limit, for q > 0.5, and move to an HD profile for 
q < 0.1. An analysis of the deviations from the fast and 
slow interconversion limits and its impact in the nonequi- 
librium phase diagram remains an interesting challenge 
which requires a systematic study. 
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FIG. 10: (Color online) Cargo density profile from Monte- 
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D± = E± = 0. (a) q —¥ oo(fast interconversion limit), (b) 
q = 2, (c) q = 0.5, (d) q = 0.25, k+ = 0.1 (e) q = 0.1 , (f) 
g = 0.05. 
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